Excitation spectra of the spin- 1/2 triangular-lattice Heisenberg antiferromagnet 
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We use series expansion methods to calculate the dispersion relation of the one-magnon excita- 
tions for the spin- 1/2 triangular- lattice nearest-neighbor Heisenberg antiferromagnet above a three- 
sublattice ordered ground state. Several striking features are observed compared to the classical 
(large-5) spin-wave spectra. Whereas, at low energies the dispersion is only weakly renormalized 
by quantum fluctuations, significant anomalies are observed at high energies. In particular, we find 
roton-like minima at special wave-vectors and strong downward renormalization in large parts of the 
Brillouin zone, leading to very flat or dispersionless modes. We present detailed comparison of our 
calculated excitation energies in the Brillouin zone with the spin- wave dispersion to order 1/5 cal- 
culated recently by Starykh, Chubukov, and Abanov cond-mat/0608002 . We find many common 
features but also some quantitative and qualitative differences. We show that at temperatures as low 
as 0.1J the thermally excited rotons make a significant contribution to the entropy. Consequently, 
unlike for the square lattice model, a non-linear sigma model description of the finite-temperature 
properties is only applicable at extremely low temperatures. 

PACS numbers: 75.10.Jm 



O 
o 



> 

oo 

o 
o 

oo 
o 

O 

-i— > 

a 



T3 

o 
o 



X 



I. INTRODUCTION 



Spin- 1/2 Heisenberg antiferromagnets on frustrated 
lattices constitute an important class of strongly corre- 
lated quantum many-body systems. The interest in these 
models has been particularly stimulated by the tantaliz- 
ing possibility that the interplay between quantum fluc- 
tuations and geometric frustration might lead to a spin- 
liquid ground state and fractionalized (i.e., S = 1/2) 
"spinon" excitations. By a spin liquid we mean a state 
which breaks neither translational nor spin rotational 
symmetry. This exotic scenario originated with the pio- 
neering work by Anderson and Fazekas more than thirty 
years ago^i where they suggested that a short-range res- 
onating valence bond (RVB) state might be the ground 
state of the nearest-neighbor Heisenberg antiferromag- 
net on the triangular lattice. More recently, consider- 
able progress has been made in understanding such states 
in terms of field theory and quantum dimer models. 2 
Whereas the existence of such a ground state and of spin- 
half excitations is well established in one dimension^ 3 , it 
is yet to be conclusively established theoretically in a re- 
alistic two-dimensional Heisenberg model^ 

Among the most important such models is the one con- 
sidered by Anderson and Fazekas, namely the Heisen- 
berg antiferromagnet on the triangular lattice with only 
nearest-neighbor (isotropic) exchange interactions (here- 
after just referred to as the triangular-lattice model 
for brevity). However, over the past decade numerical 
studies£i£iL& using a variety of different techniques do 
not support the suggestion in Ref. 1 of a spin-liquid 
ground state for this model. Instead, they provide evi- 
dence that the ground state is qualitatively similar to the 
classical one, with noncollinear magnetic Neel order with 



a three-sublattice structure in which the average direc- 
tion of neighboring spins differs by a 120 degree angle. 

On the other hand, there are other theoretical re- 
sults which suggest that some properties of this model 
are indeed quite unusual. First, a short-range RVB 
state is found to have excellent overlap with the exact 
ground state for finite systems, much better than for 
the square latticeAi^ Second, variational calculations for 
RVB states, both with and without long range order, give 
very close estimates for the ground state energies miHL 
Third, early zero-temperature series expansion studies 
found some evidence that this model may be close to 
a quantum critical point M Fourth, one can make gen- 
eral arguments, based on the relevant gauge theories^ 
that the quantum disordered phase of a non-collincar 
magnet should have deconfmed spinons^ although in 
the ordered phase the spinons are confinedii 3 " Finally, 
high-temperature series expansion studies^ performed 
for temperatures down to J/4 (J being the exchange in- 
teraction) found no evidence for the "renormalized clas- 
sical" behavior that would be expected from a semiclas- 
sical nonlinear sigma model approach, if the ground state 
has long range order ^ISilL The actual behavior (sum- 
marized in some detail in Sec. IVII Afl is rather striking 
and is in stark contrast to the square-lattice model for 
which the "renormalized classical" behavior appears very 
robust >i£ 

On the experimental side, there are currently no ma- 
terials for which it has been clearly established that their 
magnetic properties can be described by the S = 1/2 
isotropic triangular-lattice Heisenberg antiferromagnet. 
In contrast, it has been clearly established that the 
one dimensional and square lattice Heisenberg models 
with only nearest-neighbour exchange give good descrip- 
tions of a number of materials. Examples of the for- 
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mer include KCuF^ and Sr 2 Cu03r& and of the lat- 
ter Cu(DCOO) 2 .4D 2 Oi2i For the triangular-lattice the 
most exciting prospect may be the organic compound 
k-(BEDT-TTF) 2 Cu 2 (CN) 3 ^ which is estimated from 
quantum chemistry calculations to have weak spatial 
anisotropy 23 . Indeed comparisons of the thermodynamic 
susceptibility with the high temperature expansions cal- 
culated for a class of spatially anisotropic triangular- 
lattice models suggest that the system may actually be 
very close to the isotropic triangular-lattice Heisenberg 
modelpS4 In Section IVIII Bl we review the recent exper- 
imental evidence that this material has a spin liquid 
ground state. 

In this work we present series expansion calculations 
for the triangular-lattice model. The primary focus is on 
the dispersion relation of the magnon excitations above 
the 120-degree spiral-ordered ground state. A brief de- 
scription of some of our results was presented in an earlier 
communication. 25 In this paper we discuss our results 
and the series expansion and extrapolation methods in 
more detail. We also compare our results quantitatively 
with very recent calculations of the spin-wave dispersion 
by Starykh et al£& based on nonlinear spin-wave theory 
which includes quantum corrections of order 1/S (to be 
called SWT+l/S results) to the classical large- S or linear 
spin- wave theory (LSWT) results. 

One of the most striking features of the spectrum is the 
local minimum in the dispersion at the six wave vectors in 
the middle of the faces of the edge of the Brillouin zone. 
Such a minimum is absent in the spectrum calculated in 
linear spin wave theory. In particular, along the edge of 
the Brillouin zone the semi-classical dispersion is a max- 
imum, rather than a minimum at this point. This dip is 
also substantially larger than the shallow minima which 
occurs in the square lattice model. Hence, this unique 
feature seems to result from the interplay of quantum 
fluctuations and frustration. We have called this feature 
a "roton" in analogy with similar minima that occur in 
the excitation spectra of superfiuid 4 He2£ and the frac- 
tional quantum Hall effect^ In those cases, by using the 
single mode approximation for the dynamical structure 
factor one can see how the roton is associated with short 
range static correlations. Thus, an important issue is to 
ascertain whether this is also the case for the minima 
we consider here. Calculations of the static structure 
factor for the square lattice do not show a minima at 
the relevant wavevectorjS We note that the roton we 
consider is quite distinct from the "roton minima" for 
frustrated antiferromagnets that has been discussed by 
Chandra, Coleman, and Larking. The effect they discuss 
only occurs for frustrated models which have large num- 
ber of classically degenerate ground states not related by 
global spin rotations. If we apply their theory to the 
triangular lattice it does not predict such a minima. Fi- 
nally, we note that anomalous roton minima also appear 
in the spectrum of the Heisenberg model with spatially 
anisotropic exchange constants on the triangular lattice 
in the regime where the magnetic order is collincar. 25 



Such roton minima were also found in a recent study of 
an easy-plane version of the same modelj^i where the el- 
ementary excitations of the system are fermionic vortices 
in a dual field theory. In that case, the roton is a vortex- 
anti-vortex excitation making the "roton" nomenclature 
highly appropriate! 

The spectra we have calculated show substantial de- 
viations from the LSWT results, especially at high en- 
ergies and for wavevectors close to the crystallographic 
zone boundary, emphasizing the importance of non-linear 
effects in the spin dynamics. Several features of our cal- 
culated spectra are captured by the nonlinear spin-wave 
theory^ but there are also quantitative and qualitative 
differences. Both calculations show a substantial down- 
ward renormalization of the classical spectra. However, 
the highest excitation energies in the Brillouin Zone are 
lowered with respect to LSWT by about 40% in the se- 
ries calculations compared to about 25% in SWT+l/S 
results. Both calculations show substantial flat or nearly 
dispersionless spectra in large parts of the Brillouin Zone. 
However, the flat regions are much more pronounced in 
the series calculations near the highest magnon energies, 
whereas they are more pronounced at intermediate ener- 
gies in SWT+l/S results. In the series calculations the 
magnon density of states (DOS) has an extremely sharp 
peak near the highest energies, whereas in the SWT+l/S 
calculations the largest peak in DOS is at intermediate 
energies. The roton-like minima at the mid-points of 
the crystallographic zone-boundary are much more pro- 
nounced in the series calculations. They are much weaker 
in the SWT+l/S results and are really part of the flat 
energy regions contributing to the largest DOS peak in 
the latter calculations. 

The SWT+l/S calculations are much closer to scries 
expansion results than LSWT and on this basis one can 
conclude that the anomalous results obtained in series 
expansions are perturbatively related to LSWT. In other 
words, a picture based on interacting magnons captures 
the single-magnon excitations, once the non-linearities 
are taken into account. Indeed, we find that, if we treat 
magnons as non-interacting Bosons and calculate their 
entropy from the DOS obtained in the series calculations, 
we get an entropy per spin of about 0.3 at T/J = 0.3, 
a value not far from that calculated in high tempera- 
ture series expansions^ Furthermore, we find that the 
low energy magnons give the dominant contributions to 
the entropy only below T/J « 0.1. This provides a 
natural explanation for why the non-linear sigma model 
based description, which focusses only on the low energy 
magnons, must fail above T/J = 0.1. 

There remains, however, an important open question 
with regard to the spectra of this model. The question re- 
lates to the nature of the multi-particle continuum above 
the one-magnon states. In particular, how much spectral 
weight lies in the multi-particle continuum, and can it be 
described by an interacting-magnon picture, or is it bet- 
ter thought of in terms of a pair of (possibly interacting) 
spinons? This question is also related to the physical ori- 
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gin of the roton minima. We note that neutron scattering 
measurements have observed a substantial multi-particle 
continuum in the two-dimensional spin- 1/2 antiferromag- 
net Cs 2 CuCLi 3 — (which is related to the triangular lattice 
explored here, the main difference being that CS2CUCI4 
has spatially anisotropic exchange couplings). For that 
system nonlinear spin- wave theory2&24 could not account 
quantitatively for the continuum lineshapes observed ex- 
perimentally. 

The plan of the paper is as follows. In Sec. |H] we 
present our model Hamiltonian. In Sec. IIIII we dis- 
cuss the series expansion methods used for studying zero- 
temperature properties including the excitation spectra. 
Tables of various series coefficients are also presented 
there. In Sec. IIVI we discuss series extrapolation tech- 
niques. After a short Sec. about ground state proper- 
ties, we present results for the magnon dispersion in Sec. 
IV1I and compare them in detail with nonlinear spin- wave 
theory. In Sec. IVIII we consider how thermal excitation 
of the rotons affects thermodynamic properties at much 
lower temperatures than might be expected, in analogy 
with superfluid 4 He. We show how this can explain the 
absence of the renormalized classical behavior at finite 
temperatures, well below T = J/4. In Sec. IVIIII we dis- 
cuss a possible interpretation of the roton in terms of con- 
fined spinon-anti-spinon pairs and the relevance of our re- 
sults to experiments on /c-(BEDT-TTF) 2 Cu 2 (CN) 3 . Fi- 
nally, our conclusions are given in Sec. IIXI 



II. MODEL 

We consider an antiferromagnetic 5 = 1/2 Heisenberg 
model on a triangular lattice. More precisely we will 
analyze a 2-parameter Hamiltonian of this type, given 
by ' 

H = J 1 Y J S l -S n + J 2 Y,S l -S r (1) 

<in> <ij> 

Here the Si are spin- 1/2 operators. The first sum is over 
nearest-neighbor sites connected by "horizontal" bonds 
(bold lines in Fig. ^| and exchange interaction J%, the 
second sum is over nearest-neighbor sites connected by 
"diagonal" bonds (thin lines in Fig. ^) with exchange 
interaction J 2 . As far as results are concerned, in this 
paper our focus is on the isotropic model defined by J\ — 
J 2 = J, but we find it convenient to distinguish between 
J\ and J 2 for the purpose of making our discussion of the 
series expansion (Sec. Illlfl method more general. 



III. SERIES EXPANSIONS 

In order to develop series expansions for the model in 
the ordered phase, we assume that the spins order in the 
xz plane, with an angle q between neighbors along J 2 
bonds and an angle 2q along the J\ bonds. The angle q 




FIG. 1: Exchange interactions in the Heisenberg model (0 
on the triangular lattice. In this paper we focus on results for 
the case Ji — J2 = J. 



is considered as a variable; the actual value of q is that 
which minimizes the ground state energy. We rotate all 
the spins so as to have a ferromagnetic ground state, with 
the resulting Hamiltoniani 29 i 35 i 36 

H = H 1 + J X H 2 + J 2 H 3 (2) 

where 

ffi = Ji cos (2q) J2 S?Sn + J2 cos (g) ]T 5?5|, (3) 

<in> <ij> 

H 2 =J2 S?Sy+co S (2q)S*S* + sui(2q)(S*S*-S*S>), 

<in> 

(4) 

H 3 = ^ s f s j + cos (q)SfS? + sin (g)(5?5| - Sf 5|). 

<ij> 

(5) 

We introduce the Heisenberg-Ising model with Hamilto- 
nian 

H(X) =H + XV (6) 

where 

H = Hi -t^Sf -1/2), (7) 

i 

V = JiH 2 + J 2 H 3 + tJ2( S i - 1/2). (8) 

i 

The last term of strength t in both Hq and V is a lo- 
cal field term, which can be included to improve conver- 
gence. At A = 0, we have a ferromagnetic Ising model 
with two degenerate ground states. At A = 1, we arrive 
at our Heisenberg Hamiltonian of interest. We use linkcd- 
cluster methods to develop series expansion in powers of 
A for ground state properties and the magnon excitation 
spectra. The ground state properties are calculated by a 
straightforward Rayleigh-Schrodinger perturbation the- 
ory. However, the calculation of the magnon excitation 
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requires new innovations compared to a case of collinear 
order. Since S z is not a conserved quantity here due to 
the last terms in Eqs. (@J and (JjJJ, the one-magnon state 
and the ground state belong to the same sector. The 
linked-cluster expansion with the traditional similarity 
transformational fails, as it allows an excitation to anni- 
hilate from one site and reappear on another far away, vi- 
olating the assumptions for the cluster expansion to hold. 
To get a successful linked-cluster expansion, one needs 
to use the multi-block orthogonality transformation in- 
troduced in Ref. |3^ Indeed, we find that with proper 
orthogonalization the linked-cluster property holds. 

The series for ground state properties have been com- 
puted to order A 13 , and the calculations involve a list of 
4 140 438 clusters, up to 13 sites. These extend previ- 
ous calculations^ by two terms, and are given in Table 
[I] Since, we are working here with a model that has the 
full symmetry of the triangular lattice, the series for the 
magnon excitation spectra can be expressed as: 



cos(yfc ;E )cos(^y5fe 



A(fc x , k y )/J = Y2^ T ^2 Cr > r 

r— m,n 

+ cos(fc J/ v / 3(?Ti + n)/4) cos(fc2;(m — 3n)/4) 

+ cos{k y V3{m-n)/4)cos(k x (m + 3n)/4) /3 (9) 

This series has been computed to order A 9 , and the cal- 
culations involve a list of 38959 clusters, up to 10 sites. 
The series coefficients c r m n for t = 1 are given in Table 

eh 



IV. SERIES EXTRAPOLATIONS 

In this section we discuss some details of the series 
extrapolation methods used in our analysis. In order to 
get the most out of the series expansions we have adopted 
a number of strategies. The convergence of the series 
depends on the parameter t. This parameter is varied to 
find a range where there is good convergence over large 
parts of the Brillouin zone. However, the naive sum of 
the series is never accurate at points where the spectra 
should be gapless. This is true for any model and its 
reasons are explained below. 

We have found it useful to also develop series for the 
ratio of our calculated dispersion A(fc) to the classical 
(large-S) dispersion Alsw(&) obtained from linear spin- 
wave theory. Following Ref. 39, Alsw(&) for arbitrary 
A and t is given by 



A LSW (fc) = 2Sy/(\A + C){XB + C) 



(10) 



where 



A = Ji cos(fc K ) + 2J 2 cos(fca;/2) cos(v / 3fc y /2), 
B = J\ cos{k x ) cos(2g) 



,k x , , V3 
— l 

C = 2t(l - A) - Ji cos(2q) - 2J 2 cos(g). 



+ 2J 2 cos(y) cos( — k y ) cos(g), 



We can expand ALsw(k) in powers of A, and the ratio 
of our series expansion calculation A(k) to the series for 
this linear spin- wave energy ALSw(k) will be called the 
ratio series for the rest of the paper. The naive sum of 
this ratio series appears to converge better because to get 
estimates for A(k) from it, we need to multiply the sum 
by the classical energy ALsw(k) and this ensures that 
both vanish at the same k-points. 

We have also done a careful analysis of the series us- 
ing series extrapolation methods. By construction, the 
Hamiltonian H (A) has an easy-axis spin-space anisotropy 
for A < 1, which leads to a gap in the magnon disper- 
sion. This anisotropy goes away in the limit A — > 1 when 
the Hamiltonian becomes SU(2)-invariant. In this limit 
the gap must also go away as long as the ground state 
breaks SU(2) symmetry. This closing of the gap is known 
to cause singularities in the series. The singularities are 
generally weak away from ordering wavevectors and gap- 
less points, but are dominant near the ordering wavevec- 
tor where the gap typically closes in a power-law manner 
in the variable 1 — Ai 35 i 36 

We have used d-log Pade approximants and integrated 
differential approximants in our analysis. In general, 
these approximants represent the function of interest / 
in a variable x by a solution to a homogeneous or inho- 
mogeneous differential equation, usually of first or second 
order, of the form 



Pk{x) 



11 

dx 2 



) L (x)^+R M (x)f+S T (x) = (11) 
ax 



where Pk, Ql, Rm> St are polynomials of degree K, 
L, M, T respectively. The polynomials are obtained by 
matching the coefficients in the power series expansion in 
x for the above equation. They are uniquely determined 
from the known expansion coefficients of the function / 
and can be obtained by solving a set of linear equations. 
If Pk and St are set to zero, these approximants corre- 
spond to the well known d-log Pade approximants, which 
can accurately represent power-law behavior. Integrated 
differential approximants have the additional advantage 
that they can handle additive analytic or non-analytic 
terms which cause difficulties for d-log Pade approxi- 
mants. It is also possible to bias the analysis to have 
singularities at predetermined values of x with or without 
predetermined power-law exponents. Using such approx- 
imants, which enforce a certain type of predetermined 
behavior on the function, is called biased analysis. We 
refer the reader to Ref. |3g for further details. 

We found that the convergence near the ordering 
wavevector Q (see Fig. was particularly poor. We 
know that we must have gapless spectra at k = Q and 
k = as long as there is long-range order in the sys- 
tem. Yet, most unbiased analysis gave a moderate gap 
at k = Q. The convergence is better near k — 0, where 
unbiased analysis is consistent with very small values of 
the gap. This behavior near Q may be some evidence 
that long-wavelength correlations are not fully captured 
by the available number of terms in the series. 
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V. GROUND STATE PROPERTIES 

In this section we briefly discuss results for two ground- 
state properties of the triangular-lattice model: the 
ground state energy per site Eg/N and the Neel order 
parameter M (i.e., the sublattice magnetization). 

In Fig. |21 we show the extrapolated ground state en- 
ergy as a function of the angle q between nearest neig- 
bor spins along J2 bonds. Clearly the ground state 
energy is minimized when q takes the classical value 
27r/3. The resulting value for the ground state energy 
is Eo/N = —0.5502(4) J, which compares well with re- 
sults obtained from other methods (see Table ITTT)l . 

The series for the order parameter M is extrapolated 
assuming a square-root singularity at A = 1, i.e., we ex- 
trapolate the series in the variable 5 = 1 — (1 — A) 1 / 2 
using integrated differential approximants. This leads to 
the estimate M — 0.19(2). This estimate is known to be 
sensitive to the choice of the power law£ Our value for 
M shows good consistency with what is obtained from 
other methods (see Table ITTT)l . 

-0.535 I 




FIG. 2: The ground state energy per site Eo/N, as a function 
of the angle q between nearest neighbor spins along J2 bonds, 
for the triangular-lattice model (i.e., Ji = J2 = J). The 
minimum energy is obtained when q — 2tt/3, the same as for 
the classical model. 



VI. EXCITATION SPECTRA 

The triangular-lattice Brillouin zone with selected 
wavevectors is shown in Fig. [3] In Fig. 0] we plot our 
most carefully extrapolated spectra along selected direc- 
tions of the Brillouin zone using integrated differential 
approximants with appropriate biasing near the gapless 
points. The error bars are a measure of the spread in the 
extrapolated values from different approximants. Also 
shown in figure are the results from naive summation of 
the series as well as naive summation of the ratio series 
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FIG. 3: Reciprocal space of the triangular lattice including 
the hexagonal first Brillouin zone. Squares denote ordering 
wavevectors, circles denote wavevectors of the "roton" min- 
ima. The labeled points have coordinates O = (0,0), P = 
0/3,0), A = (tt,0), B = (7r,7r/V3), C = (27r/3,27r/V3), 
Q = (4tt/3, 0), and E = (0,n/y/3). Also shown is the path 
ABOCPQBE along which the magnon dispersion is plotted 
in Figs. 0] and 
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FIG. 4: (Color online) Calculated spectra along ABOCPQBE 
of the Brillouin zone. Series extrapolation results (data points 
with error bars) are plotted together with naive sum of series 
with t = 2 (green curve) and naive sum of ratio series with 
t = 1 (blue dashed) and t = 2 (magenta dashed). 



(see section IV) with two different t values. In Fig. [5] 
we plot the series expansion results together with LSWT 
and SWT+l/S spectra. 

The comparisons in Fig. 0] show that the substantial 
depression in the spin-wave energies obtained in the se- 
ries expansions, over large parts of the Brillouin zone, is 
a very robust result that does not depend on extrapola- 
tions. Roton-like minima at wavevector B is also a very 
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FIG. 7: (Color online) Projection plot showing the SWT+l/S 
magnon energies in the triangular-lattice Brillouin zone. 



FIG. 5: (Color online) Magnon spectra along ABOCPQBE 
from series expansions compared with LSWT (dashed red 
line) and SWT+l/S (green line). 



In Fig. |Sj we show the density of states (DOS) obtained 
from series analysis, LSWT and SWT+l/S. In each case 
the integrated density of states is normalized to unity. 



robust result already present in naive summation of the 
series. On the other hand, series extrapolations are es- 
sential near gapless points O, C, and Q. Even though the 
ratio series gives gapless excitations at these points, it 
does not get the spin- wave velocity right. 

Since it is tedious to perform the full analysis of the 
spectra at all points of the Brillouin zone (and not nec- 
essary for the higher energy spectra), we have instead 
carried out a more restricted D-log Pade analysis over 
the whole zone. In Fig. H3 we show a two-dimensional 
projection plot of the spectra in the full Brillouin Zone. 
The color-code is adopted to highlight the higher energy 
part of the spectra, where our results should be most re- 
liable, and minimize the variation at low energies where 
this analysis is not reliable. In Fig. the correspond- 
ing two-dimensional projection plot for the SWT+l/S 
calculations arc shown. 




FIG. 6: (Color online) Projection plot showing the magnon 
energies obtained from series expansions in the triangular- 
lattice Brillouin zone. 
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FIG. 8: (Color online) Plots of magnon density of states for 
the series expansions, LSWT and SWT+l/S spectra. 



From these plots, we make the following observations: 

1. The SWT+l/S results share many common features 
with the series expansion results. Over most of the Bril- 
louin zone the SWT+l/S results fall in between LSWT 
and series expansion results. They show that quantum 
fluctuations lead to substantial downward renormaliza- 
tion of the higher energy magnon spectra. This is in con- 
trast to unfrustrated spin models, such as square-lattice 
or linear chain models, where quantum fluctuations lead 
to increase in excitation energies. 

2. There are quantitative differences in the downward 
renormalization. The highest magnon energies are low- 
ered with respect to LSWT by about 40% in the series 
results and by about 25% in SWT+l/S results. 

3. The agreement in the low energy spectra and the 
spin-wave velocities is good when SWT+l/S results are 
compared to the biased integrated differential approxi- 
mant analysis of the series. 
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4. Both the series results and SWT+l/S results show 
relatively flat or dispersionless spectra over large parts of 
the Brillouin Zone. These lead to sharp peaks in the den- 
sity of states. However, there are some qualitative and 
quantitative differences here. In the series results the flat- 
test part of the spectra that gives rise to the largest peak 
in the DOS are near the highest energy. A second smaller 
peak in the DOS primarily gets contributions from the 
neighborhood of the roton minima. Both these regions 
are highlighted in Fig. In contrast, in SWT+l/S re- 
sults the peak in the DOS near the highest magnon en- 
ergies is much smaller. The flattest part of the spectra 
in SWT+l/S calculations are from the region near the 
roton minima. These are highlighted in Fig. 11111 



series 0.88 < EIJ <0.93 




-4 -2 2 4 



FIG. 9: (Color online) Highlight of regions in the Brillouin 
zone that contribute to the DOS peaks in the series calcula- 
tions. 



5. The roton minima at wavevector B and equiva- 
lent points are present in both series expansion results 
and in SWT+l/S results. However, they arc much more 
pronounced in the series results. A similar roton min- 
ima is seen in the square-lattice at k = (tt, 0), where 
it is also more prominent in series expansion and quan- 
tum Monte Carlo results) 44 ' 4 ^ 41 ? absent in SWT+l/S re- 
sults and barely visible in the next higher order spin- wave 
results 4£ 

6. The two-dimensional plots for both series expan- 
sions and SWT+l/S have a similar look with a central 
annular high energy region, which is separated from 6 



SWT+1/S 11 <£/J<1.2 
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SWT+1/S 0.72 < E/J <0.86 




-4 -2 2 4 



FIG. 10: (Color online) Highlight of regions in the Brillouin 
zone that contribute to the DOS peaks in the SWT+l/S cal- 
culations. 



high energy lobes by a minima in the middle. The an- 
nular region in SWT+l/S appears more circular than 
in the series results, although both have clear hexagonal 
features. The lobes also have some differences. 

7. The SWT+l/S calculations also predict finite life- 
times for spin- waves living around the center of the Bril- 
louin Zone. These have not been taken into account in 
the series calculations and may also contribute to the 
difference between the two spectra. 

Overall the comparison shows that the SWT+l/S re- 
sults have many common features but also some differ- 
ences. Neutron scattering spectra on a triangular-lattice 
material would be very exciting to compare with. In the 
meanwhile, sharp peaks in DOS may already be singled 
out in optical measurements. Given the qualitative differ- 
ences, such measurements should be able to differentiate 
bwteen the SWT+l/S and series results. However, such 
spectra would depend on various matrix elements, and 
it would be important to develop a detailed theory for 
Raman scattering for these systems. 



VII. FINITE TEMPERATURE PROPERTIES 

In this section we discuss the implications of the spec- 
tra we have calculated, particularly the rotons, to prop- 
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erties of the triangular lattice model at finite tempera- 
tures. To emphasize the importance of this issue we first 
discuss the renormalized classical behavior expected at 
low temperatures for two-dimensional quantum spin sys- 
tems, and the results of earlier high temperature series 
expansions which suggested otherwise. 

A. Finite-temperature anomalies 

For a two-dimensional quantum spin system with an 
ordered ground state, the low temperature behavior 
should correspond to a Renormalized Classical (RC) 
one, that is a "classical" state with interacting Gold- 
stone modes which is captured by the non-linear sigma 
models It is instructive to compare the behavior of 
square and triangular lattices. Spin wave theory sug- 
gests that there are not significant differences between 
the quantum corrections for square and triangular lat- 
tice models at zero temperature. For both lattices, a 
diverse range of theoretical calculations suggest that the 
reductions in sublattice magnetization and spin-stiffness 
are comparable. If this is the case then one might ex- 
pect Renormalized Classical behavior in the model to 
hold upto comparable temperatures. The relevant model 
for the square lattice is the 0(3) model, and it has been 
very successful at describing both experimental results 
and the results of numerical calculations on the lattice 
models 

For the triangular lattice, there are three Goldstone 
modes, two with velocity c± and one with velocity C|| . 
The corresponding spin stiffnesses are denoted, p± and 
p«. Several different models have been suggested to be 
relevant, including 0(3)xO(3)/0(2)tt and SU(2)i& All 
of these models predict similar temperature dependences 
for many quantities. 

In the large N expansion, including fluctuations to or- 
der 1/N (the physical model has N = 2), the static struc- 
ture factor at the ordering wavevector isi£ 

S(Q)^0.85(^-) 4 £(T) 2 (12) 

where the correlation length £(T) (in units of the lattice 
constrant) is given by^& 

f(T) = „, 21 (i)(^)' /2 « P (^) (13) 

where c = (2cj_ + C||)/3 and p s — (2p± + pn)/3 is the 
zero-temperature spin stiffness, which sets the temper- 
ature scale for the correlations. These expressions are 
quite similar to those for the 0(3) model that is relevant 
to the square lattice, with the 4ir replaced by 2ir. An 
important prediction of equations (|12fl and i|13|) is that a 
plot of Tln(S(Q)) or Tln(T£ 2 (T)) versus temperature at 
low temperatures should increase with decreasing tem- 
perature and converge to a finite non-zero value which 



is proportional to the spin stiffness in the ordered state 
at zero temperature. Indeed, the relevant plots for the 
spin- 1/2 square lattice modeli4 and the classical triangu- 
lar lattice models do show the temperature dependence 
discussed above. However, in contrast, the plots for the 
spin-1/2 model on the triangular lattice do not. In partic- 
ular, T\n(S(Q)) or Tln(T£ 2 (T)) are actually decreasing 
with decreasing temperature^ down to 0.25J. This is 
what one would expect if the ground state was actually 
quantum disordered with a finite correlation length at 
zero temperature. Hence, to be consistent with the or- 
dered ground state at zero temperature these quantities 
must show an upturn at some much lower temperature. 

The zero temperature value of the spin stiffness has 
been estimated for the TLM by a variety of methods, as 
shown in Table Ifffl The values arc in the range 0.06 J 
to 0.09J. For p s = 0.06J and c = Ja taken from non- 
linear spin wave theory (also consistent with the disper- 
sion relation found by series expansions), Equation 113|) 
implies that the correlation length should be about 0.6 
and 12 lattice constants at temperatures of T = J and 
T = 0.25 J, respectively. For comparison, the high tem- 
perature series expansions givei^ values of about 0.5 and 
1.5 lattice constants, at T = J and T = 0.25J, respec- 
tively. It should be noted that the definitions of the 
correlation length in the field theory and in the series 
expansions is slightly different 

Furthermore, the entropy for the non-linear sigma 
model at low temperatures is just that of non-interacting 
bosons in two dimensions, 



s{T) = A ( 4 + 4" J T 2 + 0{T 4 ) (14) 



where A is a dimensionless constant of 0(1). SWT+l/S 
gives C|| = 1.11J and c± = 0.69J>ii^ This means that 
for T <C J the system should have very small entropy. In- 
deed for the square lattice this is the case: it is about 0.05 
at T = 0.3J. Quantum Monte Carlo calculations found 
that for T < 0.25J the internal energy for the square 
lattice had the corresponding T 3 dependence^ How- 
ever, for the triangular lattice the entropy is still 0.3 at 
T = 0.3J. Chubukov, Sachdev, and Senthiljii suggested 
that the origin of the above discrepencies was related to 
a crossover between quantum critical and renormalised 
classical regimes. Previously, we suggested 2 ^ that the 
above discrepancies could be explained if one considered 
the rotons to be composed of a spinon and anti-spinon 
which were excited thermally. However, we now show 
how thermal excitations of rotons can explain the above 
discrepencies. There is a significant analogy here with 
the role that rotons play in superfluid He where they 
start to make substantial contributions to the entropy at 
temperatures much less than the roton gap^i^ 
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B. Contributions of rotons to finite temperature 
properties 

We will calculate the entropy of the triangular-lattice 
model by assuming that the magnon excitations can be 
treated as a gas of noninteracting bosons with a disper- 
sion as obtained from the series calculations. The entropy 
per site for noninteracting bosons (measured in units of 
ks = 1) is given by 



s(T) 



de g(e) 



e/T 

ol3e _ l 



In 



(i -.-*■) 



(15) 



where the DOS g(e) is normalized to unity. A plot of this 
entropy calculated from the series DOS in Fig. [S]is shown 
in Fig. ^] It is seen that the entropy is in fact larger 
than 0.3 at T = 0.3J and thus consistent with the high- 
temperature series data. It is also clear that the contri- 
bution to the entropy from the rotons/high-energy exci- 
tations (energy ~ 0.5 and above) starts dominating over 
the contribution from the low-energy Goldstone modes 
at a temperatures slightly above T = 0.1J which is only 
about 1/5 of the roton gap. This shows that for the tri- 
angular lattice model the presence of rotons significantly 
influences thermodynamic properties even at very low 
temperatures and thus provides an explanation of why 
the finite-temperature behavior of the triangular-lattice 
model are very different from the square lattice model. 
In particular it implies that a nonlinear sigma-model de- 
scription will be valid only at very low temperatures. 



Total Entropy 
Contribution from 0<e<0.5 
Contribution from 0.5<e<0.75 
Contribution from 0.75<e<1.0 
Contribution from 0.5<e<1.0 



r/2o. 




FIG. 11: Entropy of the triangular-lattice Heisenberg model 
due to the magnon excitations. Contributions to the en- 
tropy from different energy ranges are shown. The low energy 
magnons only dominate the entropy below T — 0.1J. 



spectra at fairly low energies (about 4-times lower than 
for the square-lattice), gives a natural explanation for 
why the results of high-temperature series expansions 
for the square and triangular lattice models are quali- 
tatively different. However, there are still some impor- 
tant questions to be addressed. First, we should stress 
that we still lack a physical picture (such as exists for 
superfluid 4 He, thanks to Feynman2!) of the nature of 
the rotons. Moreover, recent experimental results on 
k-(BEDT-TTF) 2 Cu 2 (CN) 3 , and recent variational RVB 
calculations with spinon excitations, raise further issues. 
We now briefly review these points. 



A. Are magnons bound spinon-anti-spinon pairs? 

In our earlier paper— we suggested a possible expla- 
nation of the 'roton' minima in the magnon dispersion 
relation in terms of a downward energy renormalization 
due to a level-repulsion from a higher-energy two-spinon 
continuum. This requires that the spinon dispersion have 
local minima at specific wave vectors. For the square- 
lattice model a similar interpretation, based on the 7r-fiux 
phase ^ was originally proposed by Hsu— (see also Ref. 
l57l) to explain the minima observed at (n, 0) in that case. 
For the triangular- lattice model, several RVB states have 
spinon excitations with minima at the required locations. 
For example, Lee and FengS and OgataS considered a 
Gutzwiller projected BCS state with d x 2_ y 2 + id xy pair- 
ing symmetry in the Gutzwiller approximation. How- 
ever, the variational energy of this state (— 0.484J) is 
about 15 per cent higher than the best estimates of the 
true ground state energy (see Table ITTT|) . In contrast, the 
Gutzwiller projected BCS states recently studied by by 
Yunoki and Sorella— using variational Monte Carlo have 
energies comparable to the best estimates of the ground 
state energy. In their study, a state which can be related 
to a short range RVB state, has very good variational en- 
ergy, and has spinons with mean field dispersion relation 

E(k) = (n 2 + (A) 2 (cos 2 fci + cos 2 fc 2 + sin 2 (A: 1 - fc 2 ))) 1/2 

(16) 

where k\ and fc 2 are the components of k that are parallel 
to the reciprocal lattice vectors, G\ and G 2 , of the trian- 
gular lattice. This dispersion has local minima at the four 
wave vectors, k = \ (±Gi ± G 2 ) . Spinon-antispinon ex- 
citations which make up spin triplet excitations will then 
have local minima at the six points in the middle of the 
edges of the Brillouin zone, i.e., the location of the roton 
minima found in the series expansions. 



VIII. DISCUSSION 

We have seen that SWT+l/S results for the spectra 
share many features of the series expansion results. Fur- 
thermore, the existence of rotons and flat-regions in the 



B. Experimental results 

We now review recent experimental results on the Mott 
insulating phase of k-(BEDT-TTF) 2 Cu 2 (CN) 3 . Very in- 
terestingly, this material does not show any magnetic 
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long-range order down to the lowest temperature stud- 
ied, 32 mKj2£ despite the fact that this temperature is 
four orders of magnitude smaller than the exchange in- 
teractions estimated to be around 250 K. 

The temperature dependence of the Knight shift and 
nuclear magnetic relaxation rate, 1/Ti, associated with 
13 C nuclei which have a significant interaction with the 
electron spin density have also been measured for this 
material iSLS This is a particularly useful measurement 
because the relaxation rate gives a measure of the range 
of the antiferromagnetic correlations. The observed tem- 
perature dependence of the Knight shift is the same as 
that of the uniform magnetic susceptibility^ as it should 
be. As the temperature decreases, the ratio 1/T\T in- 
creases by a factor of about two from 300 K down to 10 
K, at which it decreases by about thirty per cent down 
to 6 K. There is no sign of splitting of NMR spectral 
lines, as would be expected if long range order develops. 
In contrast, for k-(BEDT-TTF) 2 Cu 2 [N(CN) 2 ]C1 mate- 
rial, 1/TiT increases rapidly with decreasing temperature 
and exhibits a cusp at the Neel temperature, reflecting 
the diverging antiferromagnetic correlation length. Evi- 
dence for the existence of magnetic order, in the latter 
material, comes from the splitting of NMR lines at low 
temperatures^ 3 - 

The observed temperature dependence of 1/T\T and 
the spin echo rate 1/T 2 for k-(BEDT-TTF) 2 Cu 2 (CN) 3 
is distinctly different from that predicted by a non- 
linear sigma model in the renormalised classical regime 16 , 
namely that 1/TiT be proportional to T 5 / 2 £(T), and 
1/T 2 be proportional to T 3 £(T), where the correlation 
length is given by (|13fl . In particular, if this material has 
a magnetically ordered state at low temperatures, then 
both 1/T\T and 1/T 2 should be increasing rapidly with 
decreasing temperature, not decreasing. In the quan- 
tum critical regime, close to a quantum critical point 
1/Ti ~ T v where r\ is the anomalous critical exponent 
associated with the spin-spin correlation function. Gen- 
erally, for O(n) sigma models, this exponent is much less 
than one. If rj > 1, as occurs for field theories with decon- 
fined spinonsjiS^ then 1/TiT decreases with decreasing 
temperature, opposite to what occurs when the spinons 
are confined, because then 77 1 It is very interesting 
that at low temperatures, from 1 K down to 20 mK, it 
was founds that 1/Ti ~ T 3 / 2 and 1/T 2 ~ constant. In 
contrast for the materials described by the Heisenberg 
model on a square lattice 64 or a chainj2£ both relaxation 
rates diverge as the temperature decreases. Hence, these 
NMR results are clearly inconsistent with a description 
of the excitations of this material in terms of interacting 
magnons. It is also observed that a magnetic field induces 
spatially non-uniform local moments^ Motrunich has 
given a spin liquid interpretation of this observation^ 
The simplest possible explanation of why these results at 
such low temperatures are inconsistent with what one ex- 
pects for the nearest neighbour Heisenberg model is that 



such a model may not be adequate to describe this mate- 
rial and the spin liquid state may arise from the presence 
of ring exchange terms in the HamiltonianiSi This has 
led to theoretical studies of such models. 66 - 67 



IX. CONCLUSIONS 

The present comparisons of the series results of the 
spectra with the order 1/S spin- wave theory suggests 
that the shape of the magnon dispersion relation can be 
understood in a more conventional picture of interact- 
ing magnons. Furthermore, the existence of the roton 
minima at points in the middle of the edges of the Bril- 
louin zone and regions of flat dispersion in the zone, can 
explain why the low temperature properties of the tri- 
angular lattice model are so different from those of the 
square lattice. However, we still lack a clear physical pic- 
ture for the nature of the rotons. An important issue to 
resolve is whether the most natural description for them 
is in terms of bound spinon-antispinon pairs. 

From a theoretical point of view it may be interesting 
to add other destabilizing terms to the Hamiltonian, such 
as second neighbor interactions and ring-exchange terms, 
which can demonstrably lead to short correlation lengths 
and destabilize the 120 degree order, and then explore the 
changes in the dispersion relation, one-magnon weight 
and continuum lineshapes with variations in the model 
parameters. 

It is also important to examine these results in the con- 
text of the organic material k-(BEDT-TTF) 2 Cu 2 (CN) 3 . 
The measured uniform susceptibility for this material 
shows excellent parameter-free agreement with the cal- 
culated susceptibility for the spin- 1/2 triangular-lattice 
Heisenberg model. Yet, this system does not develop 
long-range order down to T/J ps 10~ 4 . 
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TABLE I: Series coefficients for the ground state energy per site Eo/N and the order parameter M for t = and t = 1 for the 
isotropic triangular-lattice model (q = 2n/3). Series coefficients of X n up to order n = 13 are listed. 
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TABLE II: Series coefficients for the magnon dispersion for the isotropic triangular-lattice model, calculated for t — 1 in Eqs. 
J7J and Nonzero coefficients c r>mjn in Eq. JUJ up to order r — 9 are listed. 
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TABLE III: Ground state properties of the isotropic triangular-lattice model, obtained by different methods. Eo/N is the 
energy per site (in units of J) for a system with N lattice sites. The order parameter M is the value of the expectation value of 
the spin in the ordered state. It would have a value of 0.5 in the absence of quantum fluctuations, and is zero in a spin liquid 
ground state. p s is the average spin stiffness (in units of J) which in a non-linear sigma model description sets the temperature 
scale of the finite temperature properties. Note that some of the spin liquid states based on variational wave functions give 
values for the ground state energy comparable to the best estimates. DMRG, QMC, V, GF, and ED denote density matrix 
renormalisation group, Quantum Monte Carlo, Variational, Greens function, and Exact diagonalization, respectively. SRVB 
denotes short range RVB. GA denotes the Gutzwiller approximation. SB+l/N denotes Schwinger boson mean-field theory 
with 1/N fluctuations. 
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